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Abstract 

This  research  models  and  analyzes  the  distribution  of  heat  and  current  in  a  buffered 
superconducting  or  hyper-conducting  wire  that  shows  potential  for  use  in  different 
capacities  in  multiple  Air  Force  systems  including  the  Active  Denial  System.  The 
thesis  includes  a  brief  background  of  the  react  ion- diffusion  system  of  partial  differen¬ 
tial  equations  provided  by  AFRL/RZPG  and  development  of  the  numerical  scheme. 
It  then  explores  solutions  to  the  model.  These  solutions  indicate  some  of  the  various 
heat-related  failures  that  may  be  observed  in  such  a  wire.  The  nature  of  the  solu¬ 
tions  observed  depends  on  the  characteristics  of  the  wire,  operating  temperature  and 
efficiency  of  cooling. 
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HEAT  AND  CURRENT  PROPAGATION  IN  BUFEERED  SUPERCONDUCTING 


AND  HYPER-CONDUCTING  WIRE 

I.  Introduction 

Over  the  last  two  decades  the  Department  of  Defense  has  shifted  its  focus  more 
and  more  towards  peace  keeping  and  addressing  asymmetric  threats.  Along  with  this 
shift  have  come  increased  hghting  in  urban  areas  and  a  need  to  minimize  collateral 
damage.  There  has  consequently  been  an  increased  need  for  non-lethal  weapons.  In 
2000  the  Air  Force  began  testing  the  Active  Denial  System,  a  less-than-lethal  crowd 
control  system  utilizing  a  direct-able  microwave  beam.  The  Air  Force  is  currently 
sponsoring  work  to  develop  an  aerial  variant  of  the  Active  Denial  System  for  use  on 
helicopters. 

Active  Denial  and  similar  systems  require  large  amounts  of  power.  In  order  to 
generate  large,  megawatt-level  amounts  of  electrical  power,  either  high  voltage  or 
large  levels  of  current  need  to  be  available.  Voltage  is  limited  at  higher  altitudes  due 
to  lower  atmospheric  pressure  [5].  Voltages  that  are  considered  low  at  sea  level  will 
produce  arcing  at  higher  altitudes.  For  that  reason,  the  way  to  increase  power  on  an 
aircraft  is  not  to  increase  voltage  but  to  increase  current.  An  increase  in  current  will 
require  larger  gauge  and/or  more  wiring.  For  small  aircraft,  two  major  concerns  are 
the  size  and  weight  of  a  prospective  system.  Power  generation  and  transport  could 
be  accomplished  with  conventional  copper  wires  on  a  C-130,  where  weight  is  not  as 
much  of  an  issue.  Conventional  power  production  and  transport  at  megawatt  levels 
would  be  too  heavy  and  bulky  to  be  used  on  a  smaller  platform,  such  as  a  helicopter. 
In  order  to  develop  a  helicopter  mounted  Active  Denial  System,  several  components 
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need  to  be  smaller  and  lighter  than  their  ground  counterparts. 

A  high-temperature  superconducting  wire  is  considerably  smaller  and  lighter  than 
a  comparable  capacity  copper  wire.  In  terms  of  weight,  even  when  cryogenics  sys¬ 
tems  are  factored  in,  high-temperature  superconductors  have  shown  potential  sav¬ 
ings  of  roughly  100  kg/m  of  cable  operating  at  10  to  15  MW  DC  [5].  This  makes 
high-temperature  superconductor-based  components  a  promising  option  for  use  in 
the  spatially  limited  environment  of  a  helicopter  or  other  aircraft. 

Reliability  is  a  concern  for  any  military  system  and  a  crucial  part  of  reliability  is 
survivability.  In  high-temperature  superconductors,  one  of  the  survivability  require¬ 
ments  is  temperature  remaining  below  a  damage  threshold.  At  this  temperature  the 
superconductor  will  be  irreversibly  damaged.  This  constitutes  a  failure  of  the  system. 
To  prevent  a  failure,  it  is  necessary  to  be  able  to  detect  as  quickly  as  possible,  by 
voltage  drop  or  temperature  variance,  any  heat  pockets  that  may  develop  and  act  to 
mitigate  them.  Figure  1  shows  a  burned  out  armature  from  a  liquid-hydrogen  cooled 
generator.  The  windings  of  the  armature  are  copper,  but  the  rotors  of  the  generator 
were  high-purity  aluminum  which  is  a  hyper-conductor.  The  burned  portions  of  the 
armature  illustrate  the  severity  of  the  damage  caused  by  such  failures. 

In  current,  state-of-the-art,  high-temperature  superconducting  wires,  the  temper¬ 
ature  may  climb  in  a  localized  region  and  not  be  detectable  until  the  system  fails 
due  to  damage  from  heat  [9].  The  Power  Generation  Branch  of  the  Power  Division, 
Propulsion  Directorate,  Air  Force  Research  Laboratory  (AFRL/RZPG)  is  research¬ 
ing  buffered  high-temperature  superconductors  and  the  impact  of  varying  interfacial 
resistivity  between  conducting  layers. 

This  research  addresses  a  general  model  for  current  and  temperature  in  a  super¬ 
conducting  wire  and  the  effect  of  increasing  interfacial  resistivity  between  conducting 
layers.  The  intent  is  to  increase  the  rate  at  which  heat  will  spread  so  that  it  is  detected 
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Figure  1.  The  armature  from  a  liquid- hydrogen  cooled  generator.  The  rotors  of  the 
generator  were  high-purity  aluminum  which  is  a  hyper-conductor.  The  windings  of 
the  armature  are  copper.  Note  the  burned  areas  visible  on  the  facing  edge  of  the 
armature  [7]. 


before  the  damage  threshold  is  reached.  While  damage  or  destruction  of  the  system 
is  the  most  dramatic  type  of  heat  related  failure,  there  are  other  less  severe  types  of 
failure  that  may  be  observed.  The  range  of  failures  from  destruction  of  the  system  to 
minor  interruption  of  current  will  be  examined.  The  hnal  possibility  addressed  is  a 
system  that  will  recover  from  a  heat  related  incident  and  maintain  stable  operation. 

Chapter  2  discusses  relevant  background  and  presents  the  model  as  a  reaction- 
diffusion  system.  Equilibrium  solutions  are  explored  and  basic  dimensional  analysis 
is  performed  on  the  system.  Chapter  3  develops  a  fast,  robust  implicit /explicit  nu¬ 
merical  scheme  for  the  model.  Chapter  4  presents  solutions  of  the  reaction-diffusion 
system  as  interfacial  resistance  and  coefficient  of  cooling  are  varied.  The  solutions  are 
qualitatively  separated  into  hve  classes  according  to  morphological  behavior.  Chapter 
5  discusses  areas  for  future  research. 
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II.  Background 


2.1  Superconductors 

Superconductivity,  discovered  in  1911,  is  a  property  displayed  by  certain  materi¬ 
als  that  provide  very  low  (or  no)  electrical  resistance  as  well  as  excluding  magnetic 
fields  at  low  temperatures.  Over  the  years  several  different  materials,  originally  met¬ 
als  but  later  ceramics  and  organic  materials,  have  exhibited  superconductivity  [12], 
Superconductors  are  currently  used  in  applications  ranging  from  magnetic  resonance 
imaging  to  particle  accelerators,  primarily  as  high-held  electromagnets. 

In  a  supercondnctor,  at  a  certain  temperatnre  called  the  critical  temperature, 
there  is  a  sndden  change  from  zero  electrical  resistivity  to  high  electrical  resistivity 
as  temperatnre  increases.  In  yttrium  barium  copper  oxide  (YBCO)  for  example,  this 
transition  takes  place  over  a  2  K  span.  Figure  2  on  page  6  [12]  illnstrates  the  change 
in  resistivity  as  a  function  of  temperature  for  YBCO.  While  there  are  several  ways  to 
classify  superconductors,  for  the  purposes  of  this  research  the  distinction  of  high  or 
low  temperature  superconductors  is  sufficient. 

YBa2Cu307_5  is  a  ceramic  that  was  one  of  the  first  materials  to  display  snpercon- 
dnctivity  at  a  temperatnre  higher  than  77  K,  the  boiling  point  of  liqnid  nitrogen.  It 
is  a  semiconductor,  but  for  certain  valnes  of  5  G  [0, 1]  it  becomes  a  snpercondnctor. 
Its  snperconducting  characteristics  vary  considerably  depending  on  concentration  of 
oxygen  with  critical  temperatnre  ranging  from  about  28  K  to  92  K.  Figure  2  shows 
a  plot  of  electrical  resistivity  as  a  function  of  temperatnre.  For  S  =  0.22,  the  critical 
temperature  is  about  77  K  and  for  5  =  0,  the  critical  temperature  jumps  to  about 
92  K.  In  addition  to  high  critical  temperature,  YBCO  has  been  shown  to  be  more 
robust  than  other  supercondnctors  with  regard  to  maintaining  snperconductivity  in 
the  presence  of  nearby  magnetic  fields,  having  been  observed  to  maintain  snpercon- 
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ductivity  at  low  temperature  in  fields  higher  than  100  T  (100  Tesla)  [15].  Even  at 
approximately  75  K  the  critical  field  of  YBCO  is  in  the  neighborhood  of  50  T  [12] . 
When  S  is  greater  than  about  0.7,  YBCO  becomes  strictly  semi-conducting  [12].  Un¬ 
less  specihcally  referring  to  S  values,  optimally-prepared  (highest  oxygen  with  5  =  0) 
YBa2Cu307_5  will  be  referred  to  as  YBCO  in  the  sequel.  Other  superconducting 
materials  have  different  critical  temperatures,  but  they  all  display  an  abrupt  change 
from  zero  resistivity  to  high  resistivity  as  the  temperature  transitions  through  the 
critical  temperature. 

Low-temperature  superconductors  operate  in  temperature  ranges  less  than  about 
23  K.  This  means  that  a  typical  cooling  method  involves  liquid  helium  which  boils 
at  roughly  4  K.  Such  a  cooling  system  is  expensive  and  bulky.  Low-temperature 
superconductors  are  also  highly  susceptible  to  heat  and  magnetic  helds.  This  limits 
the  range  of  applications  for  which  they  are  suited  and  makes  them  unsuitable  for 
most  military  applications. 

If  the  temperature  of  a  superconducting  material  rises  above  its  critical  tem¬ 
perature  or  if  a  local  magnetic  held  becomes  too  powerful,  the  material  loses  its 
superconducting  capability.  This  process  is  called  quench  and  at  that  point  the  su¬ 
perconductor  may  be  damaged  or  destroyed  unless  appropriate  protective  measures 
are  in  place.  When  a  superconductor  quenches,  a  region  called  a  normal  zone  is  cre¬ 
ated  where  there  is  greater  than  zero  electrical  resistance  [7].  The  normal  zone  will 
spread  as  the  temperature  in  that  region  rises  heating  adjacent  regions  to  quench.  In 
low-temperature  superconductors,  normal  zone  propagation  occurs  fast  enough  that 
quench  can  be  detected  and  mitigated — typically  by  shutting  down  the  system  to 
restore  temperature — before  the  temperature  reaches  the  damage  threshold.  The  low 
operating  temperature  and  susceptibility  to  magnetic  helds  makes  them,  generally, 
suitable  for  a  narrow  range  of  real  world  applications  where  conditions  can  be  carefully 
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Figure  2.  Resistivity  as  a  function  of  temperature  for  different  values  of  parameter 
S  in  YBa2Cu307_5  [12].  Notice  the  rapid  shift  from  superconducting  to  electrically 
insulating  as  temperature  increases. 


controlled  and  space  is  not  a  limiting  factor,  such  as  an  MRI  in  a  hospital. 

2.2  High-Temperature  Superconductors 

High-temperature  superconductors  maintain  zero  electrical  resistivity  up  to  ap¬ 
proximately  135  K.  This  means  that  liquid  nitrogen  which  boils  at  about  77  K  can 
be  used  to  cool  these  materials.  This  makes  high-temperature  superconductors  more 
suitable  for  many  applications.  Unfortunately,  normal  zone  propagation  is  sufficiently 
slow  in  high-temperature  superconductors  that  temperature  can  reach  the  damage 
threshold  before  there  is  any  indication  of  quench  [9].  This  makes  quench  difficult  to 
detect  and  potentially  harder  to  mitigate. 

A  superconductor  may  be  coupled  with  a  conventional  conductor  to  help  prevent 
damage/destruction  brought  on  by  quench.  In  this  case  the  electrical  current  which 
ideally  travels  through  the  superconducting  material  simply  reroutes  through  the 
conventional  conductor  until  heat  dissipation  restores  low  temperature  or  the  system 
can  be  shut  down. 

There  is  some  small  amount  of  electrical  resistance  between  the  superconducting 
and  normal  conducting  layers.  Traditionally,  it  has  been  thought  that  this  interfacial 
resistance  needs  to  be  small  [9].  AFRL/RZPG  is  currently  investigating  the  effect  of 
increasing  interfacial  resistance  on  the  speed  of  normal  zone  propagation.  If  the  speed 
of  normal  zone  propagation  could  be  made  fast  enough  to  allow  detection  and  shut¬ 
down  prior  to  system  damage,  then  high-temperature  superconductors  could  be  used 
in  a  variety  of  military  applications  including  the  electromagnets  used  for  microwave 
generation  in  the  Active  Denial  System. 
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2.3  Development  of  the  Model  from  Physical  Principles 


The  superconducting  wire  is  a  ribbon  4  mm  wide  composed  of  four  layers  as  il¬ 
lustrated  in  Figure  3.  The  hrst  layer  is  a  substrate  approximately  50  |um  thick.  This 
substrate  is  made  out  of  stainless  steel,  nickel,  or  Hastelloy  (a  high-performance  al¬ 
loy).  A  layer  of  YBCO  approximately  1  pm  thick  is  deposited  on  the  substrate  and 
a  thin  (nanometers  thick)  resistive  layer  may  be  deposited  on  top  of  this.  The  final 
layer  is  50  pm  of  copper  which  acts  as  a  stabilizer  by  providing  an  alternate  route  for 
current  to  travel  if  the  temperature  rises  above  the  critical  temperature  of  the  YBCO 
and  the  YBCO  loses  its  superconductivity  [9]. 

For  the  sake  of  a  mathematically  tractable  solution,  a  simplification  made  to  the 
mathematical  model  is  the  reduction  of  a  three-dimensional  wire  to  a  one-dimensional 
model.  This  is  reasonable  because  the  thickness  and  width  of  the  wire  are  small 
enough  that  heat  can  be  considered  approximately  uniform  in  those  directions  [9]. 
Readers  interested  in  the  complete  development  of  this  model  should  refer  to  [9]. 
The  development  in  this  thesis  focuses  on  the  elements  of  the  physical  model  that 
cannot  be  passed  over  without  loss  of  context.  The  physical  phenomenon  in  question 


Figure  3.  Cross-sectional  representation  of  the  buffered  ribbon  wire.  1  represents  the 
copper  which  acts  as  a  stabilizer  allowing  current  to  flow  when  superconductivity  is  lost. 
2  represents  the  Hastelloy  substrate  which  provides  a  base  for  depositing  the  YBCO. 
Because  of  the  relatively  high  resistance  of  the  Hastelloy  compared  to  the  copper  or 
the  YBCO,  no  current  flows  through  the  Hastelloy.  3  represents  the  thin  resistive  layer 
providing  interfacial  resistance  between  the  YBCO  and  the  copper.  4  represents  the 


YBCO. 


is  described  by  a  system  of  partial  differential  equations: 


(9  (9^ 

+  (“(^))^  +  ^))(1  “  («(a^)))^ 

+  t)  -  6»(a;,  t)o) 

^  =  “  “  /(^(a^,^))(l  -  («(a^))) 


For  readability,  this  system  will  be  written  as 


Ot  =  ^xx  +  +  /(6')(1  -  nf  +  r{u^f  -  k{9  -  Oq) 

{r{u)x)x  =  u-  f{9){l  -  u). 


(1) 

(2) 


In  either  case,  the  hrst  equation  describes  heat  and  the  second  describes  current  in 
the  wire  presented  previously.  The  variables  and  parameters  in  the  above  system 
are  listed  in  Table  1  on  page  11.  All  terms  are  dimensionless  in  this  model.  In 
practice  the  operating  temperature  9  is  greater  than  the  temperature  of  the  coolant 
6*0.  The  critical  temperature  of  the  superconductor  is  dehned  to  be  6*  =  0,  so  in 
the  absence  of  a  failure  the  operating  temperature  will  be  negative.  This  means  the 
coolant  temperature,  for  obvious  physical  reasons  lower  than  any  other  temperature 
in  the  model,  is  always  a  negative  number.  The  fraction  of  total  current  that  travels 
through  the  stabilizer  m  is  a  number  between  zero  and  one,  that  is  m  G  [0, 1],  with 
M  =  0  denoting  all  current  flows  through  the  superconductor  and  u  =  1  indicating  all 
current  flows  through  the  stabilizer. 

Following  the  development  in  [8,9],  the  hrst  terms  in  (1)  9t  =  9xx  are  the  basic 
heat  equation.  This  provides  diffusion  in  the  model.  The  following  three  terms 
+  f{9){l  —  uY  +r{uxY ,  which  are  positive,  contribute  heat.  The  term  represents 
the  heat  produced  by  current  howing  through  the  stabilizer.  Notice  when  temperature 
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is  low  enough,  there  is  no  current  traveling  through  the  stabilizer  and  this  term  will 
not  contribute  any  heat.  The  /(6*)(1  — represents  the  heat  produced  by  the  current 
in  the  superconductor.  When  temperature  is  low,  this  will  contribute  little-to-no  heat 
because  f{6)  will  evaluate  to  zero,  or  nearly  zero.  The  function  f{6)  is  detailed  in  the 
following  section.  When  all  current  is  passing  through  the  stabilizer  this  term  will  not 
contribute  heat.  The  hnal  heat  building  term  r(M2.)^  represents  the  heat  produced  by 
current  crossing  the  interface  between  the  superconductor  and  the  stabilizer.  Recall 
the  resistive  layer  between  the  superconductor  and  the  stabilizer.  This  resistance  r 
is  not  constant,  for  reasons  that  will  be  detailed  below.  The  coefficient  of  cooling 
K  is  multiplied  by  the  difference  between  operating  and  coolant  temperature  0  —  6q 
and  together  they  represent  the  heat  loss  and  cooling  forces.  The  second  equation, 
essentially  Poisson’s  equation,  derives  from  the  condition  of  charge  conservation  across 
the  vertical  thickness  of  the  wire.  The  right  hand  side,  similar  to  the  hrst  equation, 
comprises  a  term  for  current  through  the  stabilizer  u,  and  a  term  for  current  through 
the  superconductor  f{9){l  —  u). 

For  the  superconducting  wire,  6  is  dimensionless  temperature  derived  from  the 
following  relationship  (variables  are  dehned  in  Table  1): 


e 


T-Ti 

Te-Ti 


(3) 


The  current  sharing  temperature  Ti  can  be  thought  of  as  the  temperature  at  which 
current  starts  to  flow  through  the  stabilizer.  It  is  below  the  critical  temperature 
and  is  detailed  in  [9].  The  domain  is  given  by  a;  G  [— L,L],  L  is  dimensionless,  with 
either  reflecting  or  periodic  boundary  conditions.  Reflecting,  or  thermally  insulating, 
boundary  conditions  model  a  wire  of  length  2L,  while  periodic  boundary  conditions 
model  a  closed  loop  of  wire.  Mathematical  representations  of  the  boundary  conditions 
are  presented  in  the  following  chapter.  The  interfacial  resistance  along  the  ribbon  is 


10 


Table  1.  Variables  and  Terms  for  Equations  1  and  2 

Tc  critical  temperature 

Ti  current  sharing  temperature 

T  temperature 

To  ambient  (coolant)  temperature 

K  coefficient  of  cooling  (dimensionless) 

6  dimensionless  temperature 

6o  dimensionless  ambient  temperature 

f{6)  transition  function 

r  interfacial  resistance  (dimensionless) 

u  fraction  of  current  through  copper  stabilizer  (dimensionless) 

t  time  (dimensionless) 

X  length  (dimensionless) 


given  by  r{x).  Typically,  r{x)  is  a  constant  tq  along  the  length  of  the  ribbon,  bnt 
near  the  ends  it  is  benehcial  to  rednce  resistivity  for  the  pnrpose  of  simnlating  the 
wire  accurately.  In  the  wire,  current  enters  the  wire  throngh  the  copper  stabilizer. 
Cnrrent  will  need  to  overcome  the  interfacial  resistance  in  order  to  flow  into  the 
snpercondnctor.  To  make  snre  this  happens,  resistance  is  minimized  at  each  end  of 
the  wire.  To  model  this 


r{x)  =  e  +  ro  (^1  -  j 

where  e  =  0.01,  and  c  =  0.1.  This  modihcation  is  included  in  the  model  with  reflecting 
bonndary  conditions.  It  is  not  used  in  the  model  with  periodic  bonndary  conditions. 
In  that  case  interfacial  resistance  is  held  constant  r{x)  =  tq  in  the  entire  domain. 
This  allows  cnrrent  and  temperatnre  to  pass  from  L  to  —L  (throngh  the  bonndary) 
withont  interference.  In  this  way,  the  periodic  bonndary  conditions  simulate  a  length 
of  wire  with  nniform  activity  at  hxed  intervals. 
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2.4  Transition  Functions 


Of  particular  interest  in  the  study  of  this  system  is  the  function  that  models 
the  transition  of  the  superconductor  from  zero  to  high  electrical  resistance.  This 
is  represented  in  the  system  by  f{9).  A  peculiarity  of  superconducting  materials  is 
the  rapid,  or  immediate,  transition  from  zero  electrical  resistance  to  high  electrical 
resistance  at  the  critical  temperature.  Below  the  critical  temperature  the  current 
traveling  through  the  stabilizer  will  be  low  u  ~  0.  As  the  operating  temperature 
climbs,  current  through  the  stabilizer  will  increase  until  u  =  1  and  there  is  no  current 
moving  through  the  superconductor.  To  properly  model  this,  three  functions  have 
been  considered  as  the  transition  function  describing  electrical  resistivity  p  of  the 
superconductor.  The  first,  and  simplest,  is  the  scaled  Heaviside  function 

{0  lie  <l 

where  T  >  1. 

T  if  0  >  1 

As  shown  in  Figure  4,  this  is  a  step  function  which  equals  zero  below  the  critical 
temperature  and  a  constant  at  or  above  the  critical  temperature. 

The  next  function  considered  is  the  variable  ramp  shown  in  Figure  5.  Although 


Figure  4.  Scaled  Heaviside  function  for  resistivity  p  as  a  function  of  temperature  6. 
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shown  with  a  small  A,  in  practice  this  will  be  variable  depending  on  u  and  6 


0  if  6*  <  M 

<  T{e -u)/{i-u)  ifu<e<i 

r  iii<e. 


The  ramp  models  a  hnite  (2  K)  change  between  superconducting  and  insulating. 
This  function  has  been  used  to  describe  the  superconducting  relation  in  the  afore¬ 
mentioned  buffered  YBCO  wire  [8]. 

The  hnal  function  considered  is  the  exponential  function 


f{0)  = 


While  the  exponential  function  provides  a  smooth  transition  from  low  to  high  electri¬ 
cal  resistivity,  it  is  not  as  appropriate  for  the  accurate  modeling  of  the  superconducting 
material.  However,  this  model  does  describe  a  relevant  and  related  physical  device 
such  as  a  ribbon  wire  of  high  purity  aluminum  and  Hastelloy  separated  by  a  resistive 
layer  of  oxidized  aluminum.  The  Hastelloy  could  be  the  same  as  that  used  in  the 


Figure  5.  The  ramp  transition  function  for  resistivity  p  as  a  function  of  temperature 
6.  Note,  the  ramp  function  has  been  fixed  for  readability.  In  practice,  the  slope  of  the 
ramp  is  variable. 
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Figure  6.  The  exponential  transition  function  for  resistivity  p  as  a  function  of  temper¬ 
ature  6. 


buffered  YBCO  wire,  with  a  thin  (a  couple  of  nanometers  thick)  layer  of  aluminum 
deposited  and  allowed  to  oxidize  in  the  air.  Deposited  on  top  of  that  would  be  a 
layer  of  high  purity  aluminum.  Pure  aluminum  (at  least  99.99%  pure)  is  a  hyper¬ 
conductor.  This  means  that  its  electrical  resistivity  is  close  to  zero  at  sufficiently 
low  temperatures  [3].  Unlike  superconductors,  this  change  in  resistivity  is  gradual, 
ranging  from  7.55  x  10“^  pD  cm  at  20  K  to  2.733  pD  cm  at  300  K  [3].  Figure  7  on 
page  16  shows  the  relative  electrical  resistances  of  aluminum,  Hastelloy,  copper  and 
YBCO.  Notice  that  while  aluminum’s  resistance  gets  very  close  to  zero,  it  does  so 
over  a  large  temperature  range,  unlike  the  change  in  resistance  of  YBCO.  The  dimen¬ 
sionless  temperature  for  a  hyper-conducting  wire  made  of  high-purity  aluminum  and 
Hastelloy  is  slightly  different  than  that  of  a  superconducting  wire.  The  dimensionless 


temperature  6  is  given  by 
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T-Ti 

AT 


(4) 


AT  is  used  here  because  there  would  not  be  a  critical  temperature  T^.  in  the  hyper¬ 
conducting  wire,  instead  there  is  only  a  current  sharing  temperature  Ti  [7].  Unlike 
a  superconductor,  the  temperature  at  which  the  resistance  of  aluminum  becomes 
lower  than  that  of  Hastelloy  is  entirely  determined  by  their  relative  thickness.  For 
example,  consider  a  wire  such  that  the  current  sharing  temperature  Ti  =  65  K.  The 
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electrical  resistivity  of  aluminum  at  that  temperature  is  pAi  =  -1384  phi  cm  [3].  Since 
Hastelloy’s  electrical  resistivity  is  130  pfl  cm,  the  thickness  of  the  Hastelloy  would  be 
need  to  be  roughly  one  thousand  times  that  of  the  aluminum.  The  temperature  driven 
increase  in  resistance  of  aluminum,  and  the  corresponding  shift  in  current,  is  similar  to 
what  happens  in  the  superconducting  wire.  In  the  superconducting  wire  the  electrical 
resistance  of  the  YBCO  is  zero  when  T  <  and  current  will  flow  through  the  YBCO. 
When  T  >  Tc,  resistance  in  the  YBCO  will  be  much  higher  than  resistance  in  the 
copper.  At  that  point,  the  current  will  overcome  the  interfacial  resistance  between 
the  copper  and  the  YBCO,  transfer  to  the  copper  and  continue  flowing  through  the 
copper  unless  the  temperature  drops  again.  Figure  7  shows  the  change  in  resistance  as 
a  function  of  temperature  for  YBCO,  high-purity  copper,  Hastelloy  and  high-purity 
aluminum.  For  readability,  the  electrical  resistance  of  Hastelloy  is  shown  for  Hastelloy 
one  hundred  times  as  thick  as  the  aluminum.  In  this  case,  widths  of  all  materials  are 
equal  as  are  lengths.  Also,  the  resistance  of  YBCO  is  shown  at  one-twentieth  of  actual 
for  comparison  to  the  resistance  of  copper. 

Empirical  observations  suggest  the  specihc  function  is  not  critical  so  long  as  it 
provides  a  change  from  zero  or  near  zero  resistivity  to  high  resistivity  over  a  narrow 
temperature  range.  In  the  physical  wires  and  in  the  models  the  current  will  travel 
through  the  superconductor  or  the  hyper-conductor  so  long  as  the  temperature  is 
low  enough  to  allow  it.  Once  the  temperature  passes  the  critical  temperature,  the 
YBCO — or  the  hyper-conductor — becomes  a  resistor  and  the  current  will  reroute 
through  the  stabilizer  (copper  or  Hastelloy,  respectively).  For  this  reason  this  thesis 
will  address  the  case  of  the  exponential  function  as  transition  function. 
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Figure  7.  Resistance  as  a  function  of  temperature  for  Hastelloy,  YBCO,  aluminum  and 
copper  [3, 12, 14].  Resistance  of  YBCO  at  1  |j,m  thickness  is  shown  as  5%  of  the  actual 
resistance.  Aluminum  and  copper  are  50  pm  thick.  Hastelloy  must  be  100  times  thicker 
than  the  others  to  achieve  comparable  resistance.  All  materials  are  the  same  width, 
as  they  would  be  in  a  wire.  Note  that  resistivity  data  for  copper  terminates  at  50  K, 
while  resistivity  data  for  aluminum  was  available  down  to  low  temperature. 
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2.5  Constant  Eqnilibrium  Solutions 


Prior  to  numerical  modeling,  it  is  necessary  to  develop  a  fundamental  understand¬ 
ing  of  the  solutions  to  the  equilibrated  system.  Recall  the  system  of  equations  for 
temperature  6  and  current  u: 

Ot  =  Oxx  +  u‘^  +  -  uf  +  r{u^f  -  K{e  -  Oo) 

{ru^)x  =  u-  f{9){l  -  u) 


As  an  initial  exploration  of  this  system,  a  good  starting  point  is  determining  the 
stable  or  constant  equilibrium  solutions.  Consider  current  and  temperature  constant. 
Then  Ux  =  0,  9x  =  0,  9t  =  0,  and  9xx  =  0: 


+  /(^)(1  -  uf  -  k{9  -  9q) 

(5) 

u-f{9){l-u) 

(6) 

Now,  solving  (6)  for  u\ 


fjO) 

1  +  /W 


To  hud  the  equilibrium  solutions,  substitute  (7)  into  (5)  and  simplify. 


k{9-  9o) 


fjO) 

1  +  /W 


(7) 


(8) 


The  solutions  to  (8)  will  depend  on  k.  These  solutions  are  shown  in  Figure  8(a)  on 
page  19  as  well  as  Figures  8(b)  and  9(a),  where  current  m  is  a  function  of  temperature 
9.  In  those  hgures,  the  left  and  right  hand  side  of  (8)  are  plotted  representing  cooling 
forces  and  heating  forces  respectively.  The  solutions  are  indicated  at  the  intersections 
or  equilibrium  of  the  cooling  and  heating  terms.  If  temperature  9  is  greater  than  a 
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nearby  constant  equilibrium  solution,  then  cooling  will  be  stronger  than  the  heating 
forces  and  6  decreases  to  converge  to  the  equilibrium  solution.  On  the  other  hand, 
if  6  is  less  than  a  solution,  the  heating  forces  will  be  greater  than  the  cooling  and  6 
will  increase  to  converge  to  the  solution.  Therefore  these  solutions  are  stable  with 
regard  to  perturbations  in  temperature.  In  all  of  the  above  referenced  hgures,  the  low 
temperature  equilibrium  occurs  at  the  ambient  temperature  for  the  system  Oq.  In  the 
adiabatic  case,  where  k  =  0,  there  will  be  one  solution  and  it  will  be  a  stable  solution. 
Any  introduction  of  heat  (perturbation  of  the  system)  will  cause  a  shift  to  the  normal 
conducting  regime  and  temperature  will  grow  unbounded.  In  later  instances,  the 
largest  value  of  k  for  which  there  is  a  high-temperature  equilibrium  solution  will  be 
referred  to  as  k*.  When  0  <  k  <  k*,  there  will  be  two  stable  equilibrium  solutions  a 
high-temperature  and  a  low-temperature  solution.  When  k*  <  k  there  will  only  be 
a  low-temperature  solution.  In  Figure  9(b)  showing  solutions  in  the  K0-plane,  one 
can  see  that  for  k  =  0.06  there  would  be  three  intersections.  The  highest  of  these 
intersections  corresponds  to  the  highest  temperature  of  a  solution  to  this  system,  6'max 
for  the  considered  value  of  k. 

2.6  Traveling  Wave  Solution 

The  scaled  Heaviside  function  is  suited  to  one  analytic  approach  when  the  inter¬ 
facial  resistance  is  zero.  Recall  the  system  of  equations: 

Gt  =  Sxx  +  +  r{ux)‘^  +  f{9)(\  —  u)^  —  k{9  —  6q) 

{rux)x  =  u-  f{9){l  -  u) 
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(a)  Bistable  solutions 


(b)  Bistable  solutions 

Figure  8.  The  bistable  solutions,  indicated  by  the  solutions  of  (8).  For  n  >  k* ,  there  is 
only  one  stable  solution,  the  low  temperature  solution.  Note,  the  ramp  function  has 
been  fixed  for  readability.  In  practice,  the  slope  of  the  ramp  is  variable 
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K  =  0.15 


-4  -2  0  2  4 

9 

(a)  Bistable  solutions 


Figure  9.  The  bistable  solutions,  indicated  by  the  solu 
only  one  stable  solution,  the  low  temperature  solution. 


m 


1.12  0.14  0.16 

plane 


of  (8).  For  K  >  K*,  there  is 


Take  the  transition  fnnction  f{d)  =  rH{d),  where  H{d)  is  the  Heaviside  fnnction  and 
r  3>  1  is  a  scaling  parameter  and  let  r{x)  =  0.  Then  the  system  becomes: 


dt  —  Oxx  +  +  f{9){l  —  uY  —  k,{9  —  9o)  (9) 

0  =  u-f{9){l-u).  (10) 


The  second  eqnation  can  be  solved  for  u 


u  = 


f{0)  r 
i  +  f{9)  i  +  r 


H{9). 


When  K  <  Kmax,  the  solntion  has  two  stable  eqnilibria 


1  T 

9{x)  =  6*0  and  9{x)  =  6'max  =  9q  + 


K  1  +  T 


Consider  a  the  temperatnre  distribntion  which  consists  of  the  two  eqnilibrinm  so- 
Intions  separated  by  an  interface,  specihcally  6*max  on  the  left  and  9^  on  the  right  as 
shown  in  Fignre  10.  The  interface  is  not  in  eqnilibrinm  and  hence  it  will  evolve  over 
time.  A  traveling  wave  solution  merely  translates  without  changing  shape.  Consider 
the  traveling  wave  solution  9{x,  t)  =  9{x  —  ct)  =  9{^)  where  c  is  the  propagation  speed 
of  the  interface.  In  this  case 

=  9^^t  =  -c9^ 


and 


9x  —  9^^x  —  9(^  9 XX  — 


For  such  a  traveling  wave  solution  (9)  reduces  to 


-c9'  =  9”  +  +  f{9){l  -  uf  -  k{9  -  9o), 
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K 

Figure  10.  Equilibrium  solutions  0max  and  9q  connected  by  a  smooth  interface.  This 
interface  is  not  in  equilibrium  so  it  will  evolve  over  time. 

or  equivalently 

e”  +  cO'  +  u^  +  f{e){i  -  -  K{e  -eo)  =  o 

an  ordinary  differential  equation  in  This  equation  is  the  same  as 


e"  +  ce'  +  M  - 

1 

o 

=  0, 

for  6*  >  0 

(11) 

1 

+ 

1 

o 

=  0, 

for  6*  <  0 

(12) 

with  M  =  r/(l  +  r),  because  the  transition  function  f{6)  is  a  scaled  Heaviside 
function. 

The  homogeneous  problem 

6"  +  c6'  +  k6  =  0 

has  the  characteristic  equation  +  cr  —  k  =  0  which  has  solutions 

— c  +  \/c^  +  Ak  —c  —  +  Ak 

ri  =  -  and  ro  =  - . 

2  2 
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Only  Ti  will  be  considered,  since  r2  would  result  in  a  negative  speed  for  a  front  that 
is  supposed  to  move  to  the  right.  The  solutions  to  (11)  and  (12)  are 


6_{x)  =  Aq  exp{rix)  +  Bq  exp{r2x)  H - \-  Oq  for  6*  >  1  (13) 

K 

0^{x)  =  Ai  exp(ria;)  +  Bi  exp(r2a;)  +  6q  for  9  <  1.  (14) 


Where  Aq,  Ai,  Bq  and  Bi  are  arbitrary  constants.  To  ensure  continuity  and  smooth¬ 
ness  of  the  solution,  the  following  boundary  conditions  and  continuity  conditions 
apply 

1.  lim  9  =  9q 

X^OO 

2.  lim  6^  =  6'max 

X^—OO 

3.  0_(O)  =  1  =  0+(O) 

4.  0^0)  =  0^0) 

Imposing  boundary  conditions  1  and  2  on  (13)  and  (14)  forces  Bq  =  Ai  =  Q  and  the 
solution  becomes: 

M 

9_{x)  =  Aoexp(ria;)  H - \-  9q  for  6^  >  1  (15) 

K 

9^{x)  =  Bi  exp(r2a;)  +  9q  for  6*  <  1  (16) 


Applying  continuity  conditions  3  to  (15)  and  (16)  yields 


Aq  —  1 - ^0 

K 

Bi  =  1-9q 


(17) 

(18) 
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Applying  continuity  condition  4  to  the  solutions  and  solving  for  the  speed  c  gives 


c=(J-4+ _ 

V  V  M  +  Ki-i  +  eo)J  -i  +  Oo)  ■ 

The  traveling  wave  moves  to  the  right  (c  >  0)  when 

M 

^  ^  2  +  2|0o| 

The  speed  c  of  the  traveling  wave  can  be  be  approximated  by  letting  M  =  1,  and 
6q  =  —1.  Figure  11  shows  speed  of  the  traveling  wave  c  as  a  function  of  cooling  n. 
By  extension  a  plateau  solution  consisting  of  three  equilibrium  solutions  separated  by 
two  interfaces  such  as  shown  in  Figure  16(a)  on  page  40  will  grow  ii  k  <  M/2  +  2|0o| 
or  shrink  otherwise.  This  condition  on  k  appears  to  be  necessary  but  not  sufficient. 
Numerical  results  using  the  exponential  and  ramp  transition  functions  require  lower 
values  of  k  to  guarantee  a  traveling  wave  solution  in  the  case  of  no  interfacial  resis¬ 
tance. 


Figure  11.  The  speed  of  the  traveling  wave  solution  as  a  function  of  k. 
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2.7  Dimensional  Analysis 


It  is  useful  to  develop  an  understanding  of  the  importance  of  9  on  the  function 
f{9).  In  general,  whether  an  exponential  term  or  a  step  function  is  used  makes  little 
difference  in  the  types  of  solutions  observed.  The  specihcs  of  speed  and  regions  may 
change,  but  the  fundamental  solution-types — heat  gain  or  heat  loss — remain  the  same. 
It  is  possible  to  observe  the  effect  of  scaling  the  coefficient  of  6  in  the  exponential 
term  by  introducing  the  coefficient  a.  This  can  be  considered: 

9t  =  9xx  +  +  f{a9){l  -  r{ux)‘^  -  k{9  -  9o) 

{rux)x  =  u-  f{a9){l  -  u) 

Then  introduce  0  in  place  of  the  term  a9\ 

-0t  =  -(pxx  +  u^  +  /(0)(1  -  uf  +  r{Uxf  -  K{-(t)  -  -0o) 

0:0;  OL  OL 

{rUx)x  =  U-  /(0)(1  -  u) 

Now,  let  T  =  at  and  ^  =  x^/a.  The  system  of  equations  becomes: 

4>t  =  +  U^  +  /(0)(1  -  uf  +  ar{u^Y  -  -  0o) 

a{ru^)^  =  u-  /(0)(1  -  u) 

This  analysis  shows  that  changing  the  transition  function  from  f{6)  to  f{a6)  has  the 
same  effect  as  changing  interfacial  resistance  from  r  to  ar,  changing  the  coefficient 
of  cooling  (inhibitor)  from  k  to  K/a,  rescaling  time  from  t  to  at,  and  rescaling  space 
from  X  to  ^/ax,  while  keeping  f{d)  as  the  transition  function.  Hence,  it  is  sufficient 
to  examine  f{6)  over  a  range  of  values  for  k  and  r. 
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III.  Numerical  Method 


3.1  Considerations 

This  section  presents  the  development  and  details  of  the  numerical  method  used 
to  solve  the  system  of  equations 

dt  =  Oxx  +u^  +  /(6')(1  -  uf  +  -  K{d  -  Oq)  (19) 

{ru^)x  =  u- f\e){l-u)  (20) 

over  a  domain  x  G  [— T,  L\  with  either  periodic  or  reflecting  boundary  conditions. 
While  the  wire  can  be  considered  inhnitely  long  for  the  purposes  of  classical  analysis, 
numerical  analysis  requires  dehned  boundaries.  For  that  reason,  boundary  conditions 
will  be  as  follows.  For  periodic  boundary  conditions  current  and  temperature  will  be 
allowed  to  flow  across  the  boundaries 

e{t,  L)  =  e{t,  -L),  e^{t,  l)  =  -l), 

u{L)  =  u{—L),  and  Ux{L)  =  u^^—L). 

For  reflecting  boundary  conditions  current  will  be  hxed  at  the  boundaries  to  simulate 
current  injection  through  the  copper  stabilizer  and  the  boundaries  will  be  thermally 
insulating 

6x{t,—L)  =  0,  6x{t,L)  =  0,  u{—L)  =  1,  and  u{L)  =  1. 

Reflecting  boundary  conditions  are  used  to  model  a  wire  with  terminals  included. 
In  the  wire,  the  current  will  enter  through  the  copper  stabilizer  before  re-routing  to 
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the  superconductor.  Reflecting  boundary  conditions  require  the  speciflc  adaptation 


r{x)  =  e  +  ro(l-  eh°dW(2L))Vc"  j 

to  allow  for  the  change  in  resistance  necessary  for  current  injection  at  the  ends  of 
the  wire.  The  parameter  c  is  the  width  parameter  of  the  interface  between  resistance 
r{x)  ~  tq  to  r{x)  ~  e  at  the  boundaries.  Without  this  modification,  temperature 
builds  at  the  boundaries  as  a  result  of  the  current  crossing  the  interface  between  the 
superconductor  and  the  stabilizer.  Figure  12  on  the  next  page  shows  the  interfacial 
resistance  along  the  length  of  the  wire  with  reflecting  boundary  conditions.  Notice 
the  smooth  transition  from  about  tq  =  200  in  the  bulk  of  the  wire  to  almost  zero 
resistance  at  the  boundaries.  Periodic  boundary  conditions  represent  a  closed  loop  of 
wire.  In  this  case,  there  is  no  change  in  resistivity  so  r{x)  =  tq  throughout.  Using 
reflecting  or  periodic  boundary  conditions  results  in  essentially  the  same  solution  as 
long  as  the  boundaries  are  avoided. 

Initial  conditions  are  given  by 

ei0,x)  =  eo  +  ia-eo)e-^'/^'^'  (21) 

where  a  is  the  temperature  above  critical  temperature  6  =  0  and  d  is  the  width 
parameter.  Equation  (21)  produces  a  Gaussian  shape  that  models  the  distribution  of 
initial  heat. 

This  system  presents  two  difficulties  with  regard  to  a  numerical  solution.  The 
system  is  both  stiff  and  nonlinear.  The  stiffness  means  there  are  differing  time  scales. 
That  is,  the  solution  exhibits  both  rapid  changes  and  slow  changes.  An  explicit 
method  would  involve  a  Courant-Friedrichs-Lewy  (CFL)  condition.  At  <  Ax‘^/2, 
which  places  sever  restrictions  on  the  step  size  for  the  purpose  of  stability.  This 
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Figure  12.  Interfacial  resistance  over  the  length  of  [—L,L]. 

additional  restriction  would  drive  the  time  step  to  even  smaller  sizes,  further  extending 
the  time  to  compute  a  solution.  This  means  that  an  explicit  method,  which  would 
generally  be  simpler  to  implement,  would  be  impractical  [11].  To  achieve  the  desired 
level  of  accuracy,  one  would  need  to  use  very  small  step  sizes  making  computation 
slow.  An  implicit  method  such  as  a  Crank-Nicolson  permits  arbitrarily  large  step  sizes 
because  it  is  unconditionally  stable  with  respect  to  the  CFL  condition.  To  maintain 
reasonable  accuracy  At  =  0{Ax),  which  is  not  overly  restrictive. 

An  explicit  scheme  is  suitable  for  the  non-stiff  operators  in  the  equations.  An 
explicit  method  has  the  advantage  of  being  simpler  to  implement  than  an  implicit 
scheme,  so  it  should  be  faster.  The  system  is  solved  with  a  second-order  implic- 
it/explicit  (IMEX)  Adams-Bashforth/Crank-Nicolson  method.  The  implicit  portion 
takes  care  of  the  stiff  linear  terms,  while  the  explicit  portion  addresses  the  nonlinear 
term(s).  The  IMEX  Adams-Bashforth/Crank-Nicolson  method  is  absolutely  stable 
(A-stable).  This  means  that  it  converges  even  for  larger  step  sizes  while  maintaining 
second  order  accuracy  [11].  For  further  details  see  [10]  or  [11]. 
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3.2  Numerical  Model 


The  domain  is  discretized  using  a  uniform  mesh  Xj  =  jAx  for  j  =  0, 1, . . . ,  iV 
and  a  mesh  size  Ax.  An  IMEX  Adams-Bashforth/Crank-Nicolson  time  marching 
scheme  is  used  to  solve  (19)  at  each  time  step  tn  =  nAt.  Let  0”  be  the  numerical 
approximation  of  the  temperature  6{xj,tn)  and  be  the  numerical  approximation 
of  current  u{xj,tn).  The  scheme  is  given  by 


with 


and 


_  Pjn  1 

At  '  =  2 


4- 

«?=  (A A 


(22) 
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for  j  =  0,  2, . . . ,  A^.  The  variable  Oq  is  the  non-dimensionalized  ambient  temperature. 
The  current  jg  computed  by  solving  the  numerical  discretization  of  the  Poisson 

equation  (20) 


O+1/2C1+1  ~  (0+V2  +  O-1/2)  Uj  +  rj_ij2U^-i 


u,  -  /(e"+'''")(i  -  U,) 


(25) 


where  rj+1/2  =  r{xj+i/2).  When  r{x)  =  tq  (25)  simplihes  to 


^^0 


—  217j  +  Uj—i 

{Axy 
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for  j 


1,2, 


iV-1  with 


approximated  using  the  second-order  extrapolation 


p.n+l/2  _  3on  _ 

~  2^j  2^j 


The  systems  of  equations  (22)  and  (26)  are  not  closed.  Specihcally,  ghost  points 
are  necessary  on  the  domain  boundaries  {j  =  —1  and  j  =  N  +  1).  The  ghost  points 
are  eliminated  by  the  constraints  given  by  the  boundary  conditions.  To  implement 
periodic  boundary  conditions  (23)  and  (24)  j  =  0  and  j  =  N  are  replaced  by 
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where  tq  is  the  interfacial  resistance  as  previously  dehned.  To  implement  reflecting 
boundary  conditions,  the  following  constraints  are  used  to  eliminate  the  ghost  points 
at  a;_i  and  xat+i 
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and  using  the  second  order  approximation  u{x)  ~  {u{x  +  Ax)  +  u{x  —  Ax))/2 


Ui  +  U-i 
2 

Un+1  + 

2 


1 

1. 


’■o^(^  =  c'o-/(erA(i-£'o) 

=  c'"  -  /(e^'-'^Xi  -  u^). 

For  the  system  of  equations  in  which  there  is  an  exponential  transition  function, 
a  centered  difference  method  is  used  to  solve  (20)  for  u  which  is  input  into  (19). 
Because  the  equation  is  linear  with  respect  to  u,  the  operator  may  be  inverted  using 
Gaussian  elimination  and  the  result  is  used  in  the  time  marching  solution  to  (19). 
The  system  of  equations  in  which  there  is  a  ramp  as  the  transition  function  requires 
a  hxed  point  iterative  method  to  solve  (20),  because  the  ramp  function  is  nonlinear 
with  respect  to  u.  A  typical  step  size  when  exploring  general  forms  of  solutions  has 
been  Ax  =  0.25.  When  computing  specific  values,  such  as  front-speed.  Ax  =  0.125 


Table  2.  LI  error  by  step  size 


scale  1/2 

1/4 

1/8 

error  0.00548 

0.00127 

0.00025 
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has  been  used.  Table  2  on  the  preceding  page  shows  LI  error  relative  to  step  size.  LI 


error  is  computed  by 


N 


error  =  |(^(a^i)  —  0i)|Aa; 

i=l 

where  9{xi)  is  the  “exact”  solution  and  Qi  is  the  corresponding  computed  value  for  a 
given  step  size.  The  “exact”  value  is  determined  by  solving  the  system  with  a  small 
step  size  Ax  =  0.0625.  The  ratios  of  the  error  terms  conhrm  quadratic  convergence 
of  the  method. 

For  reflecting  boundary  conditions  the  second  derivative  operator  is  a  sparse  tri¬ 
diagonal  matrix.  For  periodic  boundary  conditions,  the  same  operator  is  a  sparse 
circulant  matrix.  Typical  runtime  for  a  simulation  with  L  =  200,  t  =  300,  1600 
mesh-points  in  space,  1200  mesh-points  in  time  and  reflecting  boundary  conditions 
will  be  about  22  seconds.  Typical  runtime  for  the  same  simulation  modihed  for 
periodic  boundary  conditions  is  about  67  seconds.  For  reference  purposes  these  times 
represent  performance  on  a  dual  core  2GHZ  processor. 
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IV.  Taxonomy  of  Solutions 


4.1  Map  of  Qualitative  Regions 

This  chapter  details  the  qualitative  characteristics  that  dehne  different  solutions 
to  the  system 


Ot  =  Qxx  +  u^  +  -  uf  +  r{u^f  -  K{e  -  9o) 

{r{u)^)^  =  u-  f{0){l  -u). 

The  solutions  can  be  grouped  into  seven  different  families  based  on  qualitative  dy¬ 
namics  with  six  of  these  categories  classihed  as  failures.  Recall  that  tq  is  interfacial 
resistance  and  k  is  the  cooling  coefficient.  When  a  wire  is  coiled  for  use  in  a  magnet, 
its  contact  with  the  coolant  is  limited.  This  is  different  from  a  wire  that  is  surrounded 
by  coolant.  For  these  reasons,  ro  and  k  have  been  varied  in  the  numerical  model  and 
the  results  studied.  The  initial  question  addressed  the  impact  of  interfacial  resistance 
on  speed  of  normal  zone  propagation.  The  results  of  this  research,  in  addition  to 
addressing  that  question,  show  solutions  that  are  not  so  simple. 

Figure  13(a)  on  the  next  page  shows  a  map,  in  the  rg/t-plane,  of  the  different 
regions  for  the  system  with  an  exponential  transition  function.  Figure  13(b)  shows 
a  map  for  the  system  with  a  variable  ramp  as  the  transition  function.  The  regions 
shown  in  the  maps  correspond  to  specihc  types  of  solutions.  These  maps  correspond 
to  solutions  generated  using  specihc  initial  conditions.  In  the  case  of  the  exponential 
transition  function,  the  initial  Gaussian  shape  was  centered  in  space  and  had  height 
and  width  parameters  a  =  d  =  6.  In  the  map  corresponding  to  the  ramp  transition 
function,  these  parameters  were  a  =  1.1,  d  =  Changing  these  parameters  will 
result  in  changes  to  these  maps.  The  solutions  have  been  categorized  based  on  quan- 
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Interfacial  resistance  tq 

(a)  Map  of  regions  with  exp(0)  transition  function 


Interfacial  resistance  tq 

(b)  Map  of  regions  with  ramp  transition  function 


Figure  13.  Maps  of  qualitative  regions  in  the  tqk  phase-plane.  In  figure  (a)  initial 
conditions,  see  equation  (21),  are  a  =  d  =  6;  k*  «  0.1101.  In  figure  (b)  initial  conditions 
are  a  =  1.1,  d  =  k*  =  0.495. 
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titative  as  well  as  qualitative  differences.  The  quantitative  descriptors  are  listed  in 
Section  4.1.1  . 

The  regions  can  be  described  as  follows.  The  Adiabatic  Region,  not  shown  in 
Figure  13(a),  corresponds  to  the  adiabatic  case  where  heat  increases,  but  k  =  0  so 
there  is  no  cooling.  In  this  region  6  increases  unbounded  and  the  prohle  of  6  resembles 
a  triangle  expanding  over  time.  Region  I  is  also  characterized  by  an  increase  in 
temperature,  but  there  is  an  upper  bound  for  6.  A  solution  in  this  region  resembles 
a  horizontally  expanding  plateau  similar  to  the  traveling  wave  solution  detailed  in 
section  2.6.  Solutions  in  Region  II  are  also  described  by  traveling  waves  but  as  an 
upper  bound  for  6,  they  have  small  peaks  that  form  at  the  trailing  edges  of  the 
moving  fronts.  Solutions  in  all  three  of  these  regions  display  increasing  temperature 
with  the  maximum  value  of  6  in  Regions  I  and  II  consistent  with  analysis  of  constant 
equilibrium  solutions  in  Chapter  2.  The  solutions  in  Region  III  can  be  described  as 
dissipative  solitons  [1].  Dissipative  solitons  differ  from  classical  solitons  in  that  their 
shapes  are  not  necessarily  constant,  and  their  rate  of  motion  is  not  determined  by  their 
initial  conditions.  Instead  they  can  pulsate  and  their  rate  of  motion  is  characterized 
by  system  parameters.  While  classical  solitons  may  interact,  in  dissipative  solitons 
interaction  is  limited  to  the  edges  of  their  tails  [1].  In  Region  III,  the  solitons  will 
be  generated  continuously  on  an  inhnite  domain  and  spread  along  the  length  of  the 
wire.  In  Region  IV,  the  initial  Gaussian  shape  splits  into  two  solitons  that  will 
begin  pnlsating  as  they  move  away  from  each  other  nntil  they  encounter  a  change  in 
interfacial  resistance  or  another  soliton.  At  that  point  they  will  be  repelled  and  move 
back  in  the  opposite  direction.  In  Region  V,  the  initial  Gaussian  nndergoes  a  change 
in  shape  and  size,  but  does  not  move.  The  resulting  distribntion  is  a  stationary, 
non-pnlsating  soliton.  In  Region  VI,  the  initial  Gaussian  dissipates  completely  and 
all  heat  is  absorbed  by  the  cooling  forces. 
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The  map  for  the  regions  for  solutions  of  the  system  with  ramp  transition  function  is 
qualitatively  similar  to  the  map  of  regions  for  the  system  with  exponential  transition. 
The  Adiabatic  Region  and  Regions  I  and  II  contain  solutions  that  can  be  described 
as  having  continuously  increasing  heat.  In  each  case,  the  boundary  locations  are 
determined  by  initial  conditions.  For  example,  increasing  the  width  parameter  d  in 
the  initial  Gaussian  shape  means  that  it  is  necessary  to  increase  interfacial  resistance 
in  order  to  observe  solutions  in  Region  V.  This  is  reasonable  because  in  equation  (1) 
the  gradient  of  the  current  Ux  is  multiplied  by  r  interfacial  resistance,  and  therefore 
r{uxY  can  held  at  a  hxed  value  for  appropriate  values  of  Ux  and  r. 

For  a  hxed  value  of  interfacial  resistance,  varying  the  cooling  coefficient  k  will 
change  the  types  of  solutions  that  are  observed.  For  example,  if  tq  =  20  then  k  =  0.103 
will  produce  a  solution  in  Region  I.  Increasing  cooling  so  that  k  =  0.104  produces  a 
solution  in  Region  VI.  This  qualitative  change  in  solutions  is  a  bifurcation.  While  the 
location  of  bifurcation  points  is  dependent  on  initial  conditions,  the  characteristics  of 
the  regions  are  not.  In  other  words  the  solutions  to  this  system  are  robust  with  respect 
to  background  noise  or  subtle  changes  to  initial  conditions.  6'max,  front  speed  and  the 
behavior  of  solitons  are  determined  by  the  system  conditions,  k,  r{x)  and  exp(0), 
rather  than  initial  conditions.  This  differs  from  conservative  systems  where  the  speed 
of  a  soliton  may  be  altered  by  adding  noise  or  changing  its  initial  conditions  [1]. 

As  mentioned  in  Section  2.4  the  solutions  to  the  system  of  equations  have  two 
stable  equilibrium,  normal  conductivity  corresponding  to  high  (above  critical)  tem¬ 
perature  and  superconductivity  corresponding  to  low  temperature.  In  the  first  case, 
heat  is  generated  and  there  is  not  enough  cooling  capacity  to  adequately  disperse 
it.  After  a  short  time  as  heat  builds,  this  normal  zone  will  spread  or  propagate  at  a 
constant  rate.  In  the  latter  case,  heat  is  generated,  but  it  is  quickly  dispersed.  This 
is  the  case  both  when  there  is  a  high  cooling  coefficient  k  >  k*  with  low  interfacial 
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resistance,  and  when  the  shape  of  the  initial  Gaussian  does  not  have  a  large  enough 
gradient  Ux  or  the  initial  Gaussian  is  not  high  enough  a  <  |^o|. 

Figure  14  on  the  following  page  shows  the  initial  Gaussian  shape  produced  by  (21), 
with  a  =  d  =  6,  used  to  model  initial  temperature  in  the  majority  of  this  research. 
These  are  the  initial  conditions  used  to  generate  Figure  13(a)  on  page  34.  Modifying 
the  initial  heat  distribution  can,  in  some  cases,  change  the  solution  observed.  This 
will  be  addressed  in  the  following  sections. 

4.1.1  Quantitative  Descriptors. 

Where  appropriate,  a  few  quantities  are  used  to  describe  the  solutions  and  char¬ 
acterize  the  regions.  Maximum  temperature,  referred  to  as  max(6*),  is  dehned  as 


max 

xG[—L,L] 


m)) 


and  quantihes  the  maximum  temperature  of  the  solution  in  the  domain.  Total  heat 
is  dehned  as 


9{x,  t)  dx 


'-L 


and  quantihes  the  total  temperature  above  9q.  Total  variation  is  dehned  as 


TV(t)  = 


-^9{x,  t)  \  dx. 
dx 


Total  variation  quantihes  the  total  change  in  temperature  in  the  region.  For  example, 
for  a  constant  function,  say  f{x)  =  5,  total  variation  would  be  zero.  For  f{x)  =  sin(a;) 
over  the  interval  [— 7r,7r],  total  variation  would  be  two. 
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Figure  14.  A  cross-section  of  the  initial  dimensionless  temperature  and  ratio  of  current, 
showing  the  perturbation  of  each.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this 
image. 

4.2  Adiabatic 

The  Adiabatic  case  is  the  set  of  /t  =  0  and  tq  for  which  a  temperature  profile 
6(x,t) — starting  with  the  Gaussian  distribution — will  grow  unbounded.  As  the  nor¬ 
mal  zone  propagates  along  the  length  of  the  wire,  the  temperature  will  climb  until 
current  is  shut  off  or  the  wire  is  destroyed.  Figure  15(a)  on  the  following  page  shows 
the  distribution  of  6  for  k  =  0  and  ro  =  200  at  time  t  =  62.5.  The  distribution  of  9 
resembles  a  triangle  with  height  6'max  and  base  2xq^  where  Xq  is  the  point  where  6^  =  0. 
In  the  Adiabatic  case,  total  variation  is  increasing,  as  are  total  heat  and  max(6'). 

4.3  Region  I 

Region  I  is  the  set  of  k  and  ro  for  which  a  temperature  prohle  6{x,t) — starting 
with  the  Gaussian  distribution — will  grow,  but  there  is  an  upper  bound  to  6.  The 
solution  in  this  region  resembles  a  monotonic  expanding  front,  or  a  traveling  wave.  As 
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(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  15.  (a)  A  snapshot  at  time  t  =  62.5  of  dimensionless  temperature  and  ratio  of 
current.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this  image,  (b)  The  evolution 
of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes  show  total  time  t  =  100, 
At  =  Ax  —  0.25,  red  indicates  high  temperature,  blue  indicates  low  temperature.  In 
black  and  white  images,  the  area  above/inside  each  ‘V’  is  higher  temperature  than  the 
rest  of  the  region.  In  the  adiabatic  case,  the  wave  front  is  not  as  pronounced  as  it  is  in 
Regions  I  and  II. 
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(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  16.  (a)  A  snapshot  at  time  t  =  150  of  dimensionless  temperature  and  ratio  of 
current  for  rg  =  200  and  k  =  0.08333.  The  ratio  of  current  u  G  [0, 1],  is  not  to  scale  in  this 
image,  (b)  The  evolution  of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes 
show  total  time  t  =  400,  At  =  Ax  =  0.25,  where  red  indicates  high  temperature  and  blue 
indicates  low  temperature.  In  black  and  white  images  the  area  inside  the  ‘V’  is  the 
high  temperature  region. 
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the  total  temperature  9  increases,  the  top  of  the  curve  will  flatten  out.  Figure  16(a) 
on  the  previous  page  shows  the  distribution  of  6  for  k  =  0.08333  and  tq  =  100  at  time 
t  =  150.  This  solution  corresponds  to  heat  build  up  and  the  conversion  to  normal 
mode  with  a  continuous  normal  zone  and  normal  zone  propagation.  While  there  is 
normal  zone  propagation  in  these  solutions,  the  temperature  remains  bounded  and 
should  grow  enough  to  damage  the  wire.  It  occurs  when  k  is  too  weak  to  inhibit 
quench.  For  solutions  in  this  Region,  total  variation  is  roughly  constant,  while  total 
heat  is  increasing  and  max(6*)  approaches  6'max  described  in  Chapter  2. 

4.3.1  Front  Speed. 

The  speed  of  the  traveling  wave  can  be  fonnd  nnmerically  by  tracking  the  position 
at  which  9{x,t)  =  0  and  compnting  its  speed.  Becanse  9{x,t)  =  0  lies  between 
mesh-points,  linear  interpolation  is  used  to  approximate  its  location.  Some  care  is 
taken  to  avoid  approaching  the  bonndary.  Fignre  17(a)  on  the  following  page  shows 
front-speed  as  a  fnnction  of  k  and  It  illustrates  the  fact  that  as  tq  increases,  so 
does  the  speed  of  the  traveling  wave,  as  ronghly  the  yTo,  for  valnes  of  n  as  low  as  the 
adiabatic  case.  After  compnting  front  speed,  linear  regression  is  nsed  to  determine  the 
speed  as  a  fnnction  of  k  and  rg.  In  the  system  with  exponential  transition  fnnction, 
the  speed  S  can  be  expressed  as  S{K,,ro)  =  1.26  —  16.3k  -|-  .055ro. 

While  it  is  usefnl  to  hnd  the  front-speed,  it  is  more  important  to  know  what  that 
speed  is  relative  to  the  rate  of  change  of  max(6').  Recall  that  in  the  adiabatic  or  near- 
adiabatic  case,  max(6')  may  reach  the  point  of  permanent  damage  before  a  normal 
zone  can  be  detected.  Figure  17(b)  shows  the  ratio  of  max(6*)  to  Xq  as  a  function  of 
v^- 
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(a)  Front-speed 


(b)  Ratio  of  Omax  to  Xo 


Figure  17.  (a)  Front-speed  as  a  function  of  for  parameter  0  <  k  <  0.12.  Note  that 
even  though  values  of  k  are  uniformly  spaced,  the  change  of  front-speed  is  not  uniform. 
K  =  0.12  is  the  only  value  that  is  above  the  analytically  predicted  k*  «  0.1101.  (b)  The 
ratio  of  max(0)  to  xq  as  a  function  of  for  parameter  0  <  k  <  0.006.  Note  that  k  is 
restricted  to  the  adiabatic  or  nearly  adiabatic  case. 


42 


(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  18.  (a)  A  snapshot  at  time  t  =  150  of  dimensionless  temperature  and  ratio  of 
current  for  tq  =  200  and  k  =  0.109.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this 
image,  (b)  The  evolution  of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes 
show  total  time  t  =  400,  At  =  Aa;  =  0.25,  red  indicates  high  temperature,  blue  indicates 
low  temperature.  In  black  and  white  images,  the  area  above/inside  the  ‘V’  is  the  high 
temperature  region.  In  this  region,  the  temperature  inside  the  ‘V’  is  actually  lower 
than  the  temperature  of  the  ‘V’. 
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4.4  Region  II 


Region  II  is  the  set  of  k  and  tq  for  which  a  temperature  prohle  6{x,t) — starting 
with  the  Gaussian  distribution — will  expand  as  a  traveling  wave,  in  a  manner  similar 
to  Region  I,  except  for  small  peaks  that  form  at  the  trailing  edge  of  the  moving  fronts 
of  the  solutions.  Figure  18(a)  on  the  previous  page  shows  the  distribution  of  6  for 
K  =  0.109  and  ro  =  200  at  time  t  =  150.  While  the  front  in  Region  I  is  monotonic 
in  space,  Region  II  presents  a  non-monotonic  expanding  front  with  a  local  maximum 
at  each  edge  of  the  plateau.  Temperature  9  for  the  plateau  region  is  consistent  with 
values  predicted  by  (8).  The  height  of  the  local  maximum,  6'peak,  has  been  computed 
numerically.  Figure  19  on  the  following  page  shows  the  relation  between  this  height, 
relative  to  the  predicted  ^max,  and  ro  for  several  values  of  k.  The  peaks  are  caused 
by  increased  cooling  coupled  with  high  interfacial  resistance.  The  increase  in  cooling 
makes  the  gradient  of  the  front  higher  increasing  the  heating  contributed  by  the 
term.  This  region  is  quantitatively  similar  to  Region  I;  total  variation  and  max(6') 
are  roughly  constant  with  total  heat  increasing  as  the  normal  zone  propagates.  Like 
Region  I,  the  failures  described  by  the  solutions  are  not  likely  to  lead  to  damage  or 
destruction  of  the  wire. 

4.5  Region  III 

Region  III  is  the  set  of  k  and  tq  for  which  a  temperature  prohle  6{x,t) — starting 
with  the  Gaussian  distribution — expands  outward  with  variable  6  between  the  fronts. 
The  front-speed  can  be  determined  as  for  previous  Regions.  Figure  17(a)  on  page  42 
includes  computed  speed  for  k,  and  ro  representative  of  this  Region. 

The  solution  consists  of  continuously  generated  solitons  which  will  move  unless 
they  are  kept  stationary  by  their  neighbors  and  boundary  conditions.  Since  this  is  a 
non-conservative  system  these  solitons  are  referred  to  as  dissipative  solitons  [1].  Dis- 
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Figure  19.  Overshoot  Speak  minus  predicted  maximum  Smax  as  a  function  of  for 
parameter  k  in  Region  II 

sipative  structures,  have  been  observed  in  various  biological,  chemical,  and  physical 
settings  such  as  gas-discharge  systems,  nonlinear  optical  systems,  and  nerve  pulses. 
They  have  also  been  described  by  a  variety  of  reaction-diffusion  systems  in  one  to  three 
dimensions  [1,2,4,13].  In  this  Region  normal  zones  are  created  which  subsequently 
move  but  do  not  necessarily  expand  or  dissipate.  Figure  20(a)  on  the  following  page 
shows  a  typical  prohle  of  6  for  k  =  0.115  and  tq  =  200  at  time  t  =  150. 

For  some  values  of  k,  after  the  initial  propagation  of  heat  the  solitons  will  be 
generated  until  they  have  filled  the  available  space  and  then  they  will  stop  moving. 
Figure  21(a)  on  page  48  shows  a  solution  when  k  =0.1109,  ro  =  300,  and  L  =  200. 
Notice  that  after  the  solution  reaches  the  boundaries,  there  is  no  further  change  in 
the  distribution.  In  other  cases,  solitons  will  develop  but  will  continue  pulsating  as 
shown  in  Figure  21(b),  where  k,  =0.1115,  tq  =  300,  and  L  =  200.  In  this  case,  the 
change  in  k  has  caused  the  width  of  the  solitons  to  decrease  and  they  will  continue 
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(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  20.  (a)  A  snapshot  at  time  t  =  150  of  dimensionless  temperature  and  ratio  of 
current.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this  image,  (b)  The  evolution 
of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes  show  total  time  t  =  400, 
At  =  Aa;  =  0.25,  red  indicates  high  temperature,  blue  indicates  low  temperature.  In 
black  and  white  images,  the  area  above/inside  the  ‘V’  is  the  high  temperature  region. 
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attempting  to  bud  or  generate  new  solitons.  There  is  insufficient  space  for  additional 
solitons,  so  the  solution  develops  a  pattern  in  time.  In  some  of  these  instances,  a 
pattern  of  generation  and  annihilation  develops  that  does  not  subside  over  time  as 
shown  in  Figure  22(b)  on  page  49,  where  k  =0.115,  tq  =  100,  and  L  =  210.  In 
all  cases,  solitons  will  generate  and  organize  into  a  pattern.  The  pattern  will  be 
stationary  if  length  L  is  an  integer  multiple  of  the  width  of  a  soliton.  This  pattern 
is  demonstrated  for  several  domain  sizes  in  Figures  22(a),  22(b),  22(c),  and  22(d)  on 
page  49.  As  L  increases,  the  qualitative  nature  of  the  solution,  after  reaching  the 
boundaries,  changes.  Before  the  solitons  reach  the  boundaries,  all  four  solutions  are 
the  same.  This  suggests  that  in  the  absence  of  boundaries,  in  a  very  long  wire  for 
example,  the  solitons  will  continue  moving  and  budding. 

The  nature  of  the  solutions  in  this  region  is  influenced  by  k,  tq,  and  size  of  the 
spacial  domain  L.  In  general,  the  solutions  are  similar  to  that  shown  in  Figures  21(a) 
and  21(b),  but  solutions  close  to  the  boundary  of  Region  IV  or  Region  VI — high  n — 
will  develop  fewer  solitons,  and  in  the  neighborhood  of  a  bifurcation  it  is  possible  to 
observe  symmetry  breaking  behavior. 

In  Region  III,  max(0)  is  not  predicted  by  analysis  in  Chapter  2  and  is  not  high 
enough  to  cause  damage  to  the  wire.  Solutions  in  Region  III  display  a  bounded- 
oscillating  max(6'),  increasing  total  variation  until  boundary  effects  restrict  further 
increase  and  an  increase  in  total  heat. 

4.5.1  Ramp  Transition  Function. 

Unlike  the  solutions  for  the  exponential  transition  function,  the  numerical  solu¬ 
tions  for  the  ramp  transition  function  are  sensitive  to  mesh-size  away  from  bifurcation 
points,  see  Figures  23(a),  23(b),  and  23(c)  on  page  51.  The  solutions  shown  share 
identical  values  of  k,  tq,  and  L.  The  mesh-size  decreases  from  .25  to  .0625.  In  the 
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(a)  K  =  0.1109 


(b)  K  =  0.1115 


Figure  21.  Effect  of  varying  the  cooling  coefficient  n  on  the  evolution  of  solutions  over 
time,  x-axes  show  length  when  L  =  200,  y-axes  show  total  time  t  =  6000,  At  =  Ax  =  0.25. 
(a)  A  stationary  and  non-pulsating  solution,  tq  =  300,  k  =  0.1109.  (b)  A  pulsating 

solution,  Tq  =  300,  k  =  0.1115. 
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(a)  2L  =  400  (b)  2L  =  420 


(c)  2L  =  460  (d)  2L  =  500 


Figure  22.  Effect  of  varying  the  length  of  domain  on  the  evolution  of  solutions  over 
time,  x-axes  show  length  2L  ranging  from  400  to  500,  y-axes  show  total  time  t  =  6000, 
At  =  Ax  =  0.25.  All  images  produced  with  rg  =  100,  k  =  0.115.  (a)  A  stationary,  pulsating 
system  of  solitons.  (b)  A  cycle  of  creation  and  annihilation  of  solitons.  (c)  A  stationary, 
non-pulsating  system  of  solitons.  (d)  A  stationary,  pulsating  system  of  solitons. 
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case  of  the  ramp  transition  function,  the  difference  in  computed  front-speed  due  to 
mesh-size  has  an  impact  on  the  behavior  of  the  solitons  in  this  region.  Notice  that, 
unlike  the  solutions  shown  previously,  the  difference  in  behavior  occurs  well  before 
boundary  effects  have  an  impact  on  the  solution. 

4.6  Region  IV 

Region  IV  is  the  set  of  k  and  tq  for  which  a  temperature  prohle  6{x,t) — starting 
with  the  Gaussian  distribution — will  evolve  into  traveling  solitons.  The  solitons  may 
display  pulsating  behavior  that  looks  like  a  galloping  motion  as  seen  in  Figure  25(a) 
for  parameters  k  =  0.14  and  ro  =  350.  This  region,  like  Region  III,  does  not  exist 
for  low  interfacial  resistance.  Figure  24(a)  on  page  52  shows  the  distribution  of  6  for 
K  =  0.14  and  ro  =  350  at  time  t  =  150.  This  is  similar  to  the  budding  observed  in 
Region  III,  and  the  movement  is  similar  as  well.  Each  soliton  will  start  to  bud  and 
split  into  two  solitons  but  the  cooling  force  is  strong  enough  to  prevent  a  complete 
split.  Instead,  the  soliton  will  swell,  developing  dual  peaks.  One  peak  will  die  off 
while  the  other  will  continue  in  the  direction  of  the  original  soliton.  This  budding  is 
suppressed  for  large  values  of  ro,  and  for  these  values  there  is  only  initial  movement 
to  an  equilibrium  position  as  shown  in  hgure  25(b)  where  k  =  0.19  and  ro  =  3500. 
In  the  pulsating  case,  the  solitons  move  in  one  direction  until  they  encounter  another 
soliton  at  which  point  both  will  be  repelled  and  proceed  in  the  opposite  direction. 
This  behavior  is  unlike  the  behavior  of  solitons  in  a  conservative  system.  Solitons  in  a 
conservative  system  may  interact  in  a  manner  that  includes  a  large  soliton  overtaking 
a  smaller  soliton  and  a  resulting  change  in  mass  and  speed.  In  the  case  of  dissipative 
solitons,  interaction  is  typically  limited  to  the  tails  or  edges  of  the  solitons  [1].  As  in 
Region  III,  symmetry  breaking  can  occur  in  the  bifurcation  regions.  In  considering 
the  interactions  of  these  solitons,  it  is  not  unreasonable  to  treat  them  as  particles,  or 
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(a)  At  =  Ax  =  0.25  (b)  At  =  Ax  =  0.125 


(c)  At  =  Ax  =  0.0625 


Figure  23.  Effect  of  varying  mesh-size  on  evolution  of  solutions  to  system  with  ramp 
as  transition  function,  x-axes  show  distance  L  =  150,  y-axes  show  total  time  t  =  800. 
(a)  At  =  Ax  =  0.25,  reflecting  boundary  conditions,  (b)  At  =  Ax  =  0.125,  reflecting 
boundary  conditions,  (c)  At  =  Ax  =  0.0625,  reflecting  boundary  conditions. 
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(b)  Evolution 


Figure  24.  (a)  A  snapshot  at  time  t  =  150  of  dimensionless  temperature  and  ratio  of 
current.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this  image,  (b)  The  evolution 
of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes  show  total  time  t  =  1600, 
At  =  Ax  =  0.25,  red  indicates  high  temperature,  blue  indicates  low  temperature.  In 
black  and  white  images,  the  portion  that  is  moving  in  time  is  elevated  heat.  Note, 
this  simulation  was  run  with  reflecting  boundary  conditions.  In  the  case  of  periodic 
boundary  conditions,  the  solitons  are  still  repelled  at  the  boundary,  but  they  are  closer 
to  the  edges  when  they  are  repelled. 
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(a)  ro  =  350 


(b)  ro  =  3500 


Figure  25.  Effect  of  varying  interfacial  resistance  ,  x-8Lxes  show  distance  L  =  200,  y- 
axes  show  total  time  t  =  6000,  At  =  Ax  =  0.25.  (a)  A  pulsating  soliton  solution,  tq  =  350, 
K  =  0.14,  reflecting  boundary  conditions,  (b)  A  non-pulsating  soliton  solution,  ro  =  3500, 
K  =  0.19,  reflecting  boundary  conditions. 
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as  particle-like  entities  [2,13].  Observing  the  repulsion  of  two  or  the  more  complex 
patterns  that  arise  as  multiple  solitons  move  and  are  moved  around,  is  very  similar 
to  the  behavior  of  particles  colliding,  complete  with  a  suggestion  of  conservation  of 
momentum.  Solutions  in  this  region  display  oscillating  and  bounded  values  of  max(6'), 
total  heat,  and  total  variation. 

4.7  Region  V 

Region  V  is  the  set  of  k,  and  tq  for  which  a  temperature  profile  d{x,t) — starting 
with  the  Gaussian  distribution — will  reshape  and  then  resolve  into  a  stationary  soli- 
ton.  This  reshaping  begins  with  the  peak  of  the  soliton  dropping  faster  than  the 
sides.  This  produces  a  distribution  of  6  with  two  peaks  that  will  draw  together.  Fig¬ 
ure  26(a)  on  the  following  page  shows  the  distribution  of  6  for  k  =  0.17  and  ro  =  400 
at  time  t  =  150  after  the  stationary  soliton  has  formed.  Before  the  peaks  touch,  the 
midpoint  of  the  soliton  will  push  back  up  and  stabilize.  The  heating  and  cooling 
forces  are  balanced  in  such  a  way  that  a  stable  stationary  soliton  is  sustained  in  the 
original  position.  The  interfacial  resistance  coupled  with  the  steep  gradient  term  in 
the  r{uxY  term  of  (1)  is  strong  enough  that  it  prevents  the  otherwise  powerful  cooling 
term  due  to  relatively  large  k  from  overcoming  the  weak  heating  terms.  The  solutions 
in  this  region  display  constant  values  of  max(0),  total  heat  and  total  variation  after 
the  initial  evolution. 

4.8  Region  VI 

Region  VI  is  the  set  of  k  and  tq  for  which  a  temperature  prohle  6{x,t) — starting 
with  the  Gaussian  distribution — decays  to  a  uniform  distribution  given  by  the  ambient 
temperature  6{x,t)  =  Oq.  Physically,  this  means  that  a  perturbation  in  temperature 
vanishes.  The  region  is  characterized  by  large  cooling  parameters  k  and  moderate 
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(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  26.  (a)  A  snapshot  at  time  t  =  150  of  dimensionless  temperature  and  ratio  of 
current.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this  image,  (b)  The  evolution 
of  the  solution  in  time,  x-axes  show  distance  L  =  200,  y-axes  show  total  time  t  =  400, 
At  =  Aa;  =  0.25,  red  indicates  high  temperature,  blue  indicates  low  temperature.  In 
black  and  white  images,  the  vertical  bar  in  the  middle  of  the  image  is  elevated  heat. 
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(a)  Current  and  Temperature  Distribution 


(b)  Evolution 


Figure  27.  (a)  A  snapshot  at  time  t  =  42.5  of  dimensionless  temperature  and  ratio  of 
current.  The  ratio  of  current  u  G  [0, 1]  is  not  to  scale  in  this  image,  (b)  The  evolution  of 
the  solution  in  time,  x-axes  show  distance  =  400  (x/It),  y-axes  show  total  time  =  100 
(7t),  At  =  Ax  =  0.25,  red  indicates  high  temperature,  blue  indicates  low  temperature. 
In  black  and  white  images,  the  portion  that  disappears  over  time  is  elevated  heat. 
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to  low  interfacial  resistivity  tq.  This  is  the  only  region  that  does  not  show  a  failure 
due  to  heating,  because  given  sufficient  cooling  capacity,  any  heat  that  develops  is 
diffused  and  absorbed  by  the  coolant.  Figure  27(a)  on  the  previous  page  shows  the 
temperature  for  k  =  0.132  and  r  =  200  at  time  t  =  42.5  after  it  has  dropped  below 
the  current  sharing  temperature  6  =  0.  Because  k  is  large,  the  cooling  term  K{d  —  Oq) 
dominates  the  heating  terms  of  (1)  near  the  center  of  the  Gaussian  distribution  where 
the  temperature  is  highest.  This  causes  the  peak  to  decay  quickly.  Because  the 
interfacial  resistance  is  sufficiently  strong,  the  higher  current  gradients  on  either  sides 
of  the  peak  contributes  to  heating  which  counteracts  the  cooling,  resulting  in  the 
dual  peak  distribution  shown  in  Figure  27.  The  temperature  continues  to  decay  in 
time  asymptotically  approaching  a  uniform  distribution  at  the  ambient  temperature. 
For  a  lower  interfacial  resistivity,  k  =  0.18  and  tq  =  50  for  example,  the  heat  diffuses 
more  uniformly  in  space  and  there  is  only  a  subtle  change  in  shape  of  the  temperature 
prohle  from  the  Gaussian  distribution  as  it  shrinks  until  it  disappears. 
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V.  Conclusion  and  Areas  for  Further  Research 


5.1  Conclusion 

This  thesis  has  explored  the  similarities  and  differences  between  a  model  for  heat 
and  current  in  a  superconducting  wire  and  a  similar  model  for  a  hyper-conducting 
wire.  In  each  case  there  are  solutions  characterized  by  heat  gain  and  heat  loss. 
The  solutions  characterized  by  heat  gain  describe  a  type  of  failure  in  the  wire.  By 
increasing  interfacial  resistance  between  the  superconductor  and  the  stabilizer  it  is 
possible  to  increase  the  likelihood  of  detecting  one  of  these  failures  before  critical 
damage  is  done  to  the  system.  The  differences  in  the  models  considered  are  largely  a 
matter  of  scale  in  terms  of  both  spacial  domain  and  time.  In  the  superconducting  wire, 
the  transition  to  normal  conductivity  occurs  over  a  much  smaller  temperature  range 
than  in  the  hyper-conducting  wire.  This  changes  the  scale  and  temperature  ranges  of 
some  of  the  solutions  but  does  not  change  the  qualitative  nature  of  the  solutions.  All 
research  conducted  in  conjunction  with  this  thesis  was  limited  to  the  one  dimensional 
model  with  constant  interfacial  resistance.  Expanding  these  parameters  provides 
material  for  further  research. 

5.2  Varying  Interfacial  Resistance 

It  is  possible  to  vary  the  interfacial  resistance  in  the  domain  x  G  [—L,  L].  Primarily 
this  system  has  been  studied  with  a  uniform  resistivity  over  the  domain  with  lowest 
resistivity  at  the  boundaries.  Instead,  consider  the  interfacial  resistance  as  a  series  of 
wells  in  an  otherwise  high  level  of  interfacial  resistance.  This  is  a  reasonable  simulation 
of  a  change  to  the  physical  wire.  In  the  original  case,  consider  a  continuous  wire  with 
uniform  resistive  layer  applied  between  the  YBCO  and  the  copper.  This  case  would 
be  analogous  to  placing  several  breaks  in  the  resistive  layer  along  the  length  of  the 
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wire.  Such  a  change  would  likely  slow  the  speed  of  normal  zone  propagation,  but  it 
might  also  allow  for  better  dissipation  and  more  stability  since  the  resistivity  is  one 
of  the  factors  that  contributes  to  normal  zone  propagation.  The  original  expression 
for  interfacial  resistance  is 

r{x)  =e  +  ro(l-  ehodW(2L))Vc"  j  _ 

The  simplest  approach  to  modify  r{x)  is  to  change  the  coefficient  of  x  in  the  cosine 
term 

r{x)  =  e  +  ro  (l  -  j 

where  n  G  Z’*'. 

5.3  Two-Dimensional  Model 

The  one  dimensional  model  is  sufficient  for  initial  exploration  and  describing  a 
length  of  wire.  It  is  not  appropriate  for  modeling  a  coil  of  wire  like  the  ones  used 
for  electromagnets.  In  that  case,  since  there  is  heat  transfer  between  layers,  but  no 
current  between  layers,  the  model  would  need  to  be  modified  and  or  extended. 

5.4  Transition  Function 

Consider  the  resistance  of  aluminum  shown  in  Figure  7  on  page  16.  Rather  than 
exponential,  this  appears  asymptotically  linear.  With  that  in  mind,  a  hyperbolic 
function  might  be  more  appropriate.  It  might  prove  illuminating  to  implement  a 
function  like 

/(0)  =  _  1 

instead  of  the  exponential  function. 
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